#!/usr/bin/perl -w
use strict;
use FindBin;
use Getopt::Long;
use lib "/net/cpp-group/Leo/bin";
use parse_fasta;
use Cwd;
use run_cmd;
use die_error;

#888888888888888888888888888888888888888888888888888888888888888888888888888888

#   Usage

#888888888888888888888888888888888888888888888888888888888888888888888888888888
my $usage = <<"USAGE";

USAGE:

  codeml_parallel.pl 

      DESCRIPTION:
              Runs site specific codeml analyses in parallel on the cluster. 
  
      REQUIRED PARAMETERS:
              --cdna CDNA.FILE       N.B. These should be the full paths to the files
              --prot PEPTIDE.FILE          
  
      OPTIONAL PARAMETERS:
              --f61                  Use full codon table instead of F3x4
			        
                                          N.B. This is more biologically accurate 
                                          but may over-parameterise for short sequences
              --job_name  JOBNAME    Specify this to be able to 
                                          view job progress / kill jobs later
              --remove_gaps          Remove all columns which contain ANY
                                          missing residues (i.e. a gap)
              --msf                  Use alternative formats for peptide
              --clustal                   multiple alignment 
                                          (defaults to FASTA)
              --regex "(.*)"         Regular expression to find 
                                          corresponding peptide and cDNA
                                          identifiers. (See following.)
              --verbose              Print more information
              --working_dir          Perform all analysis here instead of current directory
              --help                 This screen
  

NOTES:
    Run this programme on the cluster head node (currently fgu205 / fgu200).
        (i.e. run "ssh fgu205" to logon etc.)
  

    The corresponding cDNA and peptide identifiers must match.
        By default, identifiers are the first letters until a space/tab.
        Otherwise, the identifiers can be selected using a regular 
        expression (in --regex)

CODEML
    Nested models 1 and 2 is a simplistic model.
    Nested models 7 and 8 uses a beta distribution and 
        may be over parameterised.
    Both should predict +ve selection to be trusted.

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    Make sure positively selected sites have a omega >> 1
       or they may merely be instances of neutral evolution.

    Always check for errors in the pipeline (All files ending in '.err')
       ***ESPECIALLY***  mapping.err
	   

RESULTS:
    See results.parameters

FAILURES:
    If it fails, it may be because the tree is too large/small.
        Check the pairwise dS values in 'codeml.pairwise/pairwise.dS'.
        The corresponding fitch-predicted dS tree is in 'fitch.tree'.

   Happy Hunting.
                  .... Leo (Goodstadt) 1/11/2004 to 17/08/2007

USAGE

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#   Get options

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
# mandatory


# optional parameters
my $help = undef;
my $remove_gaps = undef;
my $regex = '(\\S+)';

my $job_name = '';
my $working_dir;
my $prot_file;
my $cdna_file;
my $msf;
my $clustal;
my $f61;
my $verbose;
my $vverbose;
{
    open CMDLINE, ">cmd.line" or die "Error:\n\tCould not open cmd.line\n";
    print CMDLINE join (" ", @ARGV), "\n";
    close CMDLINE;
}
GetOptions(
			'prot=s'	=> \$prot_file,
			'cdna=s'	=> \$cdna_file,
			'regex=s'	=> \$regex,
			'msf'		=> \$msf,
			'clustal'	=> \$clustal,
			'f61'		=> \$f61,
			'verbose'	=> \$verbose,
			'vverbose'	=> \$vverbose,
			'help'		=> \$help,
			'remove_gaps'=> \$remove_gaps,
			'job_name=s'=> \$job_name,
			'working_dir=s'=> \$working_dir,
			);

die $usage if ($help);
die $usage unless ($prot_file && $cdna_file);

my $bin_dir = $FindBin::Bin;
$working_dir ||= getcwd();

die "Error\n\t:[$working_dir/$prot_file] does not exist.\n" unless (-e $prot_file);
die "Error\n\t:[$working_dir/$cdna_file] does not exist.\n" unless (-e $cdna_file);
my $verbose_cmd = $vverbose ? "--verbose" :"";
unless ($regex =~ /^".*"$/)
{
	$regex = '"'.$regex.'"';
}

$f61 = defined ($f61) ? 0 : 3;

$verbose = 1 if ($vverbose);














#_________________________________________________________________________________________
# 
# 	Safe make directory
# 
#_________________________________________________________________________________________
sub my_mkdir($)
{
	my ($dir_name) = @_;
	if (!-e $dir_name)
	{
		mkdir $dir_name
				or die "Error\n\tCould not create the directory [$dir_name]\n$!\n";
	}
	elsif (! -d  "$dir_name")
	{
		die "Error\n\t[$dir_name] is a file not a directory.\n$!\n";
	}

}


sub print_pairwise_ctl;

print "\n\n" if $verbose;


my_mkdir ($working_dir);


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888


#		Main logic


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888


#_________________________________________________________________________________________
# 
# 	Prepare sequences
# 
#_________________________________________________________________________________________
#
#	Prepare FASTA format protein sequence file
#

sub prepare_fasta_files
{
	if ($msf)
	{
		print
	#STDERR
	 "MSF --> FASTA ...\n" if $verbose;
	run_cmd(<<"CMD");
		$bin_dir/msf_to_fasta.pl $prot_file
			> $working_dir/$prot_file.fa
CMD
	}
	elsif ($clustal)
	{
		print
	#STDERR
	 "CLUSTAL --> FASTA ...\n" if $verbose;
	run_cmd(<<"CMD");
		$bin_dir/clustal_to_fasta.pl $prot_file
			> $working_dir/$prot_file.fa
CMD
	}
	else
	{
	run_cmd(<<"CMD");
		cp $prot_file $working_dir/$prot_file.fa
CMD
	}
	run_cmd(<<"CMD");
		cp $cdna_file $working_dir/$cdna_file.fa
CMD
}




#_________________________________________________________________________________________
# 
# 	Map cDNA onto aa
# 		
# 		Can run locally
# 
#_________________________________________________________________________________________
sub map_prot_to_cdna
{
	print "Using protein alignment to align corresponding cDNA ...\n" if $verbose;
	run_cmd(<<"CMD");
		$bin_dir/map_cdna_onto_aa
			--peptide_sequences $working_dir/$prot_file.fa
			--cdna_sequences $working_dir/$cdna_file.fa
			--peptide_output $working_dir/$prot_file.mapped.fa
			--cdna_output $working_dir/$cdna_file.mapped.fa
			--regex $regex
			--missing_identifiers $working_dir/missing.cdna_identifiers.err
			--err_log $working_dir/mapping.err
			--x_for_aa_mismatches X
			$verbose_cmd
CMD
		unlink ("$prot_file.fa");
	
	if (-f "$working_dir/mapping.err")
	{
		open_or_die(*ERROR, "$working_dir/mapping.err");
		print $_ for (<ERROR>);
	}
	if (-f "$working_dir/missing.cdna_identifiers.err")
	{
		open_or_die(*ERROR, "$working_dir/missing.cdna_identifiers.err");
		print $_ for (<ERROR>);
	}
	
	#
	#	Prepare cDNA onto aa
	#
	{
		print
	#STDERR
	 "Convert mapped cDNA to PAML format ...\n" if $verbose;
	run_cmd(<<"CMD");
		$bin_dir/seq_extract
			--input_file $working_dir/$cdna_file.mapped.fa
			--wrap_line 0
			--output_format p
			$verbose_cmd
			--e $working_dir/extract.err
			> $working_dir/seq.tmp
CMD
	}
	
	if (-f "$working_dir/extract.err")
	{
		open_or_die(*ERROR, "$working_dir/extract.err");
		print $_ for (<ERROR>);
	}
	
}


my @start_val = ("0.1", "0.3", "1.3");
sub calculate_pairwise_ds
{
	#_________________________________________________________________________________________
	# 
	#
	#  -------Create pairwise Ks directory
	#
	#
	#  -------make directories
	#
	my_mkdir("$working_dir/codeml.pairwise");
	for my $start_val(@start_val)
	{
		print "Make directories for codeml analyses with a starting omega of $start_val...\n" if $verbose;
		my_mkdir ("$working_dir/codeml.$start_val");
	}
	
	
	
	
	
	
	#
	#  -------Create PAML control files ..
	#
	print "Create PAML control files ...\n" if $verbose;
	# pre-declaration
	sub print_site_specific_ctl($);
	for my $start_val(@start_val)
	{
		print_site_specific_ctl($start_val);
	}
	print_pairwise_ctl();
	
	
	
	
	#_________________________________________________________________________________________
	
	#_________________________________________________________________________________________
	#
	#  -------Calculate pairwise Ks to build tree
	#
	
	
	print "Calculate pairwise Ks to build tree ...\n" if $verbose;
	run_queue_cmd(<<"CMD", $job_name."_seq_pairs", "short_jobs.q");
	$bin_dir/seq_pairs_kaks
		-aA
		-t $working_dir/codeml.pairwise
		-c $working_dir/$cdna_file.mapped.fa
		--err_log $working_dir/kaks.pairwise.err
		-p $working_dir/paml.err
		$verbose_cmd
		> $working_dir/kaks.results
CMD
	
	if (-f "$working_dir/kaks.pairwise.err")
	{
		open_or_die(*ERROR, "$working_dir/kaks.pairwise.err");
		print $_ for (<ERROR>);
	}

}


sub build_fitch_tree_from_pairwise_ds
{
	print "Filter out Ks for build tree ...\n" if $verbose;
	run_cmd(<<"CMD");
		$bin_dir/seq_pairs_filter_kaks_results
			-S
			-R $working_dir/filter.summary
			-r 0
			-f 0
			-s 0
			-t 0
			-n 0
			-l 0
			-e $working_dir/filter.err
			$verbose_cmd
			< $working_dir/kaks.results
			> $working_dir/codeml.pairwise/pairwise.dS
CMD

	if (-f "$working_dir/filter.err")
	{
		open_or_die(*ERROR, "$working_dir/filter.err");
		print $_ for (<ERROR>);
	}
	


	#
	#  -------Build fitch tree from ks values
	#
	print "Build fitch tree from ks values ...\n" if $verbose;

	# run fitch
	run_queue_cmd(<<"CMD", $job_name."_tree_fitch", "short_jobs.q");
	$bin_dir/tree_fitch -i p
			$verbose_cmd
			< $working_dir/codeml.pairwise/pairwise.dS
			> $working_dir/fitch.tree
CMD
}






#
#  -------site specific analyses
#
sub site_specific_analyses
{
	my @run_cmds;
	for my $i (0, 1, 2)
	{
		my $description = $verbose ? "Site specific analyses with a starting omega of $start_val[$i] ...\n": undef;
		my $dir = "$working_dir/codeml.$start_val[$i]";
	
		my $codeml_path = "$bin_dir/../../tools/paml/codeml ";
		my $ctl_file = "site_specific_kaks.$start_val[$i].ctl ";
		my $output_file = "> ../output.$start_val[$i]";
	
		my $exec_cmd = $codeml_path . $ctl_file . $output_file;
		push(@run_cmds,  [$exec_cmd, $job_name . "_" . $start_val[$i], $dir, $description]);
	}

	run_queue_parallel(\@run_cmds, "medium_jobs.q", 7);

	print
#STDERR
 "Analyse results ...\n" if $verbose;

	run_cmd(<<"CMD");
	$bin_dir/parse_codeml_results.pl
			--para $working_dir/results.parameters
			--positive $working_dir/results.positive.fa
			$verbose_cmd
			$working_dir/site_specific.results.*
CMD

}







#
#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#	parameter files
#







sub print_pairwise_ctl
{
	open PAIR_CTL, ">$working_dir/codeml.pairwise/pairwise_kaks.ctl"
			or die "Error\n\t:Could not write to pairwise_kaks.ctl\n$!\n";
	print PAIR_CTL <<"CMD";
seqfile = ../seq.tmp
outfile = results
icode = 0	  			# universal code
noisy = 0
verbose = 0
runmode = -2  			# pairwise
seqtype = 1	  			# codons
CodonFreq = $f61		# F3X4; Could be 3 for codon table Was error fixed thanks to Andreas
aaDist = 0	  			# ??? equal Could be geometric or linear or G1974 Miyata
model = 0	  			# one dN/dS over tree
NSsites = 0	  			# one dN/dS over sequence
Mgene = 0		    	# ??? rates; 1: separate
fix_omega = 0
omega = .4
fix_alpha = 1		 	# gamma distribution shape set to 1 for single rate
alpha = 0		     	# gamma distribution shape set to 1 for single rate
fix_rho = 1
rho = 0
Malpha = 0		     	# ??? different alphas for genes
ncatG = 1		     	# # of categories in dG of NSsites models5D
clock = 0		     	# ??? no clock 1:clock
getSE = 1
RateAncestor = 1	 	# ??? rates alpha > 0 or ancestral rates 1 or 2
Small_Diff = .5e-6
method = 0		     	# simultaneous, 1: one branch at a time
fix_blength = 0		 	# ignore branch lengths from supplied ree
fix_kappa = 0
kappa = 2
CMD
	if ($remove_gaps)
	{
		print PAIR_CTL <<"CMD";
cleandata = 1
CMD
	}
	else
	{
		print PAIR_CTL <<"CMD";
cleandata = 0
CMD
	}

}

sub print_site_specific_ctl($)
{
	my ($start_omega) = @_;
	my $file = ">$working_dir/codeml.$start_omega/site_specific_kaks.$start_omega.ctl";
	open FITCH_CTL, $file
			or die "Error\n\t:Could not write to $file\n$!\n";
print FITCH_CTL <<"CMD";
seqfile = ../seq.tmp
treefile = ../fitch.tree
outfile = ../site_specific.results.$start_omega
icode = 0				* universal code
noisy = 0
verbose = 1
runmode = 0				* user supplied tree
seqtype = 1				* codons
CodonFreq = $f61		* 2 = F3x4,	3 = F64
aaDist = 0				* ??? equal Could be geometric or linear or G1974 Miyata
model = 0				* one dN/dS over tree
NSsites = 1 2 7 8
						* 0:one w;1:neutral;2:selection; 3:discrete;4:freqs;
						* 5:gamma;6:2gamma;7:beta;8:beta&w;9:beta&gamma;
						* 10:beta&gamma+1; 11:beta&normal>1; 12:0&2normal>1;
						* 13:3normal>0
						* one dN/dS over sequence
Mgene = 0				* ??? rates; 1: separate
fix_omega = 0
omega = $start_omega
fix_alpha = 1			* gamma distribution shape set to 1 for single rate
alpha = 0				* gamma distribution shape set to 1 for single rate
fix_rho = 1
rho = 0
Malpha = 0				* ??? different alphas for genes
ncatG = 10				* * of categories in dG of NSsites models5D
clock = 0				* ??? no clock 1:clock
getSE = 1
RateAncestor = 0		* ??? rates alpha > 0 or ancestral rates 1 or 2
Small_Diff = .5e-6
method = 0				* simultaneous, 1: one branch at a time
fix_blength = 0			* ignore branch lengths from supplied ree
fix_kappa = 0
kappa = 2CMD
CMD
	if ($remove_gaps)
	{
		print FITCH_CTL <<"CMD";
cleandata = 1
CMD
	}
	else
	{
		print FITCH_CTL <<"CMD";
cleandata = 0
CMD
	}
}


sub print_fitch_parameters
{
	my $file = ">$working_dir/fitch.parameters";
	open FITCH_CTL, $file
			or die "Error\n\t:Could not write to $file\n$!\n";
	print FITCH_CTL <<"CMD";
build_tree.dS
0
7
L
J
33
50
2
Y
CMD
	close FITCH_CTL;
}










#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888


#		main


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888


prepare_fasta_files();
map_prot_to_cdna();
calculate_pairwise_ds();
print_fitch_parameters();
build_fitch_tree_from_pairwise_ds();
site_specific_analyses();


print "All finished...\n\n";
print "Results in results.parameters\n\n";
print "Nested models 1 and 2 is a simplistic model.\n";
print "Nested models 7 and 8 uses a beta distribution and may be over parameterised.\n";
print "I normally require all 3 to predict +ve selection before I believe anything.\n";
print "If it fails, it may be because the tree is too large/small.\n";
print "   Check the pairwise dS values in 'codeml.pairwise/pairwise.dS'.\n";
print "   The corresponding fitch-predicted dS tree is in 'fitch.tree'.\n";
print "Make sure any identified +vely selected sites have a large omega >> 1\n";
print "   or they may merely be instances of neutral evolution.\n";
print "Always check for errors in the pipeline (All files ending in '.err')\n";
print "   ***ESPECIALLY*** \n";
print "   mapping.err\n\n";
print "Happy Hunting.\n\n";
print "          .... Leo (Goodstadt) 27/05/2005\n\n\n";









